## 可视化计算生态位重叠结果:
## 以0.5度buffer为例：
### 准备数据：
library(ecospat)
library(ade4)
library(raster)
library(rworldmap)
## 加载物种分布数据：#####
as_occ <- read.csv("D:/XH/xh/as.csv")[,2:3]
sa_occ <- read.csv("D:/XH/xh/sa.csv")[,2:3]
## 加载环境背景数据：
# 加载bg采样点：
#### 方法1： bg_enm D:\XH\fifth_bg2\bg_buffers
as_e_bg <- read.csv("D:/XH/fifth_bg2/bg_buffers/bg_as0.5.csv")[,2:3]
sa_e_bg <- read.csv("D:/XH/fifth_bg2/bg_buffers/bg_sa0.5.csv")[,2:3]
## 加载栅格背景数据：
env <- stack(list.files("D:/XH/third_env/envsV5tif",pattern = "tif",full.names = T))
### 构建数据组：
## 南美与亚洲：
env_as_occ <- data.frame(raster::extract(env,as_occ))
env_sa_occ <- data.frame(raster::extract(env,sa_occ))
## 方法1：bg_enm：
env_as_e <- data.frame(raster::extract(env,as_e_bg))
env_sa_e <- data.frame(raster::extract(env,sa_e_bg))
## 构建nat和inv：
## 构建nat_sa：
spocc <- c(rep(1,nrow(sa_occ)),rep(0,nrow(env_sa_e )))
envs_enm <- rbind(env_sa_occ,env_sa_e)
names(sa_e_bg) <- names(sa_occ)
sa_bg_occ <- rbind(sa_occ,sa_e_bg)
nat_sa <- na.exclude(cbind(sa_bg_occ,envs_enm ,spocc))
## 构建inv_as:
spocc3 <- c(rep(1,nrow(as_occ)),rep(0,nrow(env_as_e )))
envs_enm3 <- rbind(env_as_occ,env_as_e)
names(as_e_bg) <- names(as_occ)
as_bg_occ3 <- rbind(as_occ,as_e_bg)
inv_as <- na.exclude(cbind(as_bg_occ3,envs_enm3 ,spocc3))
names(inv_as) <- names(nat_sa)
### sa_as_enmbg_PCA ####
pca.env <-dudi.pca(rbind(nat_sa,inv_as)[3:13], center = T, scale = T, scannf = F, nf = 2)
### COMP_NICEH ####
## predict the scores on the PCA axes
scores.globclim <- pca.env$li
# PCA scores for the species native distribution
scores.sp.nat <- suprow(pca.env,nat_sa[which(nat_sa[,14]==1),3:13])$li
# PCA scores for the species invasive distribution
scores.sp.inv <- suprow(pca.env,inv_as[which(inv_as[,14]==1),3:13])$li
# PCA scores for the whole native study area
scores.clim.nat <- suprow(pca.env,nat_sa[,3:13])$li
# PCA scores for the whole invaded study area
scores.clim.inv <- suprow(pca.env,inv_as[,3:13])$li
# calculation of occurence density
z1<- ecospat.grid.clim.dyn(glob=scores.globclim,
glob1=scores.clim.nat,
sp=scores.sp.nat, R=100,
th.sp=0.05)
z2 <- ecospat.grid.clim.dyn(glob=scores.globclim,
glob1=scores.clim.inv,
sp=scores.sp.inv, R=100,
th.sp=0.05)
library(ecospat)
plot(z1$z.uncor)
plot(z1$z.cor)
ecospat.plot.niche(z1,cor = FALSE)
ecospat.plot.niche(z2,cor = FALSE)
## 绘制可视化重叠：####
ecospat.plot.niche.dyn(z1, z2, quant=0.25, interest=2,
title= "SA vs AS Niche Overlap", name.axis1="PC1(tem)",name.axis2="PC2(pre)")
ecospat.shift.centroids(scores.sp.nat, scores.sp.inv, scores.clim.nat, scores.clim.inv)
## 优化配色方案：
# col1 = colorRampPalette(c("lightpink", "indianred"),bias=1,alpha=0.4)(n=20)
# col2 = colorRampPalette(c("lightyellow", "orange"),bias=1)(n=20)
source("C:/Users/admin/Desktop/ENM_VIEW.R")
dyn(z1, z2, quant=0.25, interest=1,
title= "SA vs au Niche Overlap", name.axis1="",name.axis2="",
colZ1= "#0000004c",colZ2="#0000004c")
## 注意上述的dyn的函数结果：
dyn <- function (z1, z2, quant, title = "", name.axis1 = "Axis 1",
name.axis2 = "Axis 2", interest = 1, colz1 = "#00FF0050",
colz2 = "#FF000050", colinter = "#0000FF50",
colZ1 = "green3", colZ2 = "red3")
{
if (is.null(z1$y)) {
R <- length(z1$x)
x <- z1$x
xx <- sort(rep(1:length(x), 2))
y1 <- z1$z.uncor/max(z1$z.uncor)
Y1 <- z1$Z/max(z1$Z)
if (quant > 0) {
Y1.quant <- quantile(z1$Z[which(z1$Z > 0)], probs = seq(0,
1, quant))[2]/max(z1$Z)
}
else {
Y1.quant <- 0
}
Y1.quant <- Y1 - Y1.quant
Y1.quant[Y1.quant < 0] <- 0
yy1 <- sort(rep(1:length(y1), 2))[-c(1:2, length(y1) *
2)]
YY1 <- sort(rep(1:length(Y1), 2))[-c(1:2, length(Y1) *
2)]
y2 <- z2$z.uncor/max(z2$z.uncor)
Y2 <- z2$Z/max(z2$Z)
if (quant > 0) {
Y2.quant <- quantile(z2$Z[which(z2$Z > 0)], probs = seq(0,
1, quant))[2]/max(z2$Z)
}
else {
Y2.quant = 0
}
Y2.quant <- Y2 - Y2.quant
Y2.quant[Y2.quant < 0] <- 0
yy2 <- sort(rep(1:length(y2), 2))[-c(1:2, length(y2) *
2)]
YY2 <- sort(rep(1:length(Y2), 2))[-c(1:2, length(Y2) *
2)]
plot(x, y1, type = "n", xlab = name.axis1, ylab = "density of occurrence")
polygon(x[xx], c(0, y1[yy1], 0, 0), col = colz1, border = 0)
polygon(x[xx], c(0, y2[yy2], 0, 0), col = colz2, border = 0)
polygon(x[xx], c(0, apply(cbind(y2[yy2], y1[yy1]), 1,
min, na.exclude = TRUE), 0, 0), col = colinter, border = 0)
lines(x[xx], c(0, Y2.quant[YY2], 0, 0), col = colZ2,
lty = "dashed")
lines(x[xx], c(0, Y1.quant[YY1], 0, 0), col = colZ1,
lty = "dashed")
lines(x[xx], c(0, Y2[YY2], 0, 0), col = colZ2)
lines(x[xx], c(0, Y1[YY1], 0, 0), col = colZ1)
segments(x0 = 0, y0 = 0, x1 = max(x[xx]), y1 = 0, col = "white")
segments(x0 = 0, y0 = 0, x1 = 0, y1 = 1, col = "white")
seg.cat <- function(inter, cat, col.unf, col.exp, col.stab) {
if (inter[3] == 0) {
my.col = 0
}
if (inter[3] == 1) {
my.col = col.unf
}
if (inter[3] == 2) {
my.col = col.stab
}
if (inter[3] == -1) {
my.col = col.exp
}
segments(x0 = inter[1], y0 = -0.01, y1 = -0.01, x1 = inter[2],
col = my.col, lwd = 4, lty = 2)
}
cat <- ecospat.niche.dyn.index(z1, z2, intersection = quant)$dyn
inter <- cbind(z1$x[-length(z1$x)], z1$x[-1], cat[-1])
apply(inter, 1, seg.cat, col.unf = "#00FF0050",
col.exp = "#FF000050", col.stab = "#0000FF50")
}
if (!is.null(z1$y)) {
z <- t(as.matrix(z1$w + 2 * z2$w))[, nrow(as.matrix(z1$z.uncor)):1]
z1$Z <- t(as.matrix(z1$Z))[, nrow(as.matrix(z1$Z)):1]
z2$Z <- t(as.matrix(z2$Z))[, nrow(as.matrix(z2$Z)):1]
if (interest == 1) {
col1 = colorRampPalette(c("white", "black"),bias=1,alpha=0.8)(n=100)
col2 = colorRampPalette(c("white", "black"),bias=1)(n=100)
image(x = z2$x, y = z2$y, z = t(as.matrix(z2$z.uncor))[,
nrow(as.matrix(z2$z.uncor)):1], col = col2,
zlim = c(1e-05, cellStats(z2$z.uncor, "max")))
image(x = z1$x, y = z1$y, z = t(as.matrix(z1$z.uncor))[,
nrow(as.matrix(z1$z.uncor)):1], col = col1,
zlim = c(1e-05, cellStats(z1$z.uncor, "max")),
xlab = name.axis1, ylab = name.axis2, add = TRUE)
}
title(title)
contour(x = z1$x, y = z1$y, z1$Z, add = TRUE, levels = quantile(z1$Z[z1$Z >
0], c(0, quant)), drawlabels = FALSE, lty = c(1,
2), col = colZ1,lwd=3)
contour(x = z2$x, y = z2$y, z2$Z, add = TRUE, levels = quantile(z2$Z[z2$Z >
0], c(0, quant)), drawlabels = FALSE, lty = c(1,
2), col = colZ2,lwd=3)
contour(x = z2$x, y = z2$y, z, add = TRUE,col="gray40",lwd=1.5,drawlabels = FALSE)
}
}
arow <- function (sp1, sp2, clim1, clim2, col = "##FFFFFF00")
{
if (ncol(as.matrix(sp1)) == 2) {
arrows(median(sp1[, 1]), median(sp1[, 2]), median(sp2[,
1]), median(sp2[, 2]), col = "gray27", lwd = 4,
length = 0.1)
arrows(median(clim1[, 1]), median(clim1[, 2]), median(clim2[,
1]), median(clim2[, 2]), lty = "11", col = col,
lwd = 2, length = 0.1)
}
else {
arrows(median(sp1), 0.025, median(sp2), 0.025, col = "black",
lwd = 4, length = 0.1)
arrows(median(clim1), -0.025, median(clim2), -0.025,
lty = "11", col = col, lwd = 2, length = 0.1)
}
}
ecospat.plot.niche2 <- function (z, title = "", name.axis1 = "Axis 1", name.axis2 = "Axis 2",
cor = FALSE)
{
if (is.null(z$y)) {
R <- length(z$x)
x <- z$x
xx <- sort(rep(1:length(x), 2))
if (cor == FALSE)
y1 <- z$z.uncor/max(z$z.uncor)
if (cor == TRUE)
y1 <- z$z.cor/max(z$z.cor)
Y1 <- z$Z/max(z$Z)
yy1 <- sort(rep(1:length(y1), 2))[-c(1:2, length(y1) *
2)]
YY1 <- sort(rep(1:length(Y1), 2))[-c(1:2, length(Y1) *
2)]
plot(x, y1, type = "n", xlab = name.axis1, ylab = "density of occurrence")
polygon(x[xx], c(0, y1[yy1], 0, 0), col = "grey")
lines(x[xx], c(0, Y1[YY1], 0, 0))
}
if (!is.null(z$y)) {
col1 = colorRampPalette(c("white", "black"),bias=1,alpha=0.8)(n=100)
col2 = colorRampPalette(c("white", "black"),bias=1,alpha=0.1)(n=100)
if (cor == FALSE)
image(x = z$x, y = z$y, z = t(as.matrix(z$z.uncor))[,
nrow(as.matrix(z$z.uncor)):1], col = col2,
zlim = c(1e-06, cellStats(z$z.uncor, "max")),
xlab = name.axis1, ylab = name.axis2)
if (cor == TRUE)
image(x = z$x, y = z$y, z = t(as.matrix(z$z.uncor))[,
nrow(as.matrix(z$z.uncor)):1], col = gray(100:0/100),
zlim = c(1e-06, cellStats(z$z.cor, "max")),
xlab = name.axis1, ylab = name.axis2)
z$Z <- t(as.matrix(z$Z))[, nrow(as.matrix(z$Z)):1]
contour(x = z$x, y = z$y, z$Z, add = TRUE, levels = quantile(z$Z[z$Z >
0], c(0, 0.5)), drawlabels = FALSE, lty = c(1, 2),col = "#33FF99b2")
}
title(title)
}
dyn <- function (z1, z2, quant, title = "", name.axis1 = "Axis 1",
name.axis2 = "Axis 2", interest = 1, colz1 = "#00FF0050",
colz2 = "#FF000050", colinter = "#0000FF50",
colZ1 = "green3", colZ2 = "red3")
{
if (is.null(z1$y)) {
R <- length(z1$x)
x <- z1$x
xx <- sort(rep(1:length(x), 2))
y1 <- z1$z.uncor/max(z1$z.uncor)
Y1 <- z1$Z/max(z1$Z)
if (quant > 0) {
Y1.quant <- quantile(z1$Z[which(z1$Z > 0)], probs = seq(0,
1, quant))[2]/max(z1$Z)
}
else {
Y1.quant <- 0
}
Y1.quant <- Y1 - Y1.quant
Y1.quant[Y1.quant < 0] <- 0
yy1 <- sort(rep(1:length(y1), 2))[-c(1:2, length(y1) *
2)]
YY1 <- sort(rep(1:length(Y1), 2))[-c(1:2, length(Y1) *
2)]
y2 <- z2$z.uncor/max(z2$z.uncor)
Y2 <- z2$Z/max(z2$Z)
if (quant > 0) {
Y2.quant <- quantile(z2$Z[which(z2$Z > 0)], probs = seq(0,
1, quant))[2]/max(z2$Z)
}
else {
Y2.quant = 0
}
Y2.quant <- Y2 - Y2.quant
Y2.quant[Y2.quant < 0] <- 0
yy2 <- sort(rep(1:length(y2), 2))[-c(1:2, length(y2) *
2)]
YY2 <- sort(rep(1:length(Y2), 2))[-c(1:2, length(Y2) *
2)]
plot(x, y1, type = "n", xlab = name.axis1, ylab = "density of occurrence")
polygon(x[xx], c(0, y1[yy1], 0, 0), col = colz1, border = 0)
polygon(x[xx], c(0, y2[yy2], 0, 0), col = colz2, border = 0)
polygon(x[xx], c(0, apply(cbind(y2[yy2], y1[yy1]), 1,
min, na.exclude = TRUE), 0, 0), col = colinter, border = 0)
lines(x[xx], c(0, Y2.quant[YY2], 0, 0), col = colZ2,
lty = "dashed")
lines(x[xx], c(0, Y1.quant[YY1], 0, 0), col = colZ1,
lty = "dashed")
lines(x[xx], c(0, Y2[YY2], 0, 0), col = colZ2)
lines(x[xx], c(0, Y1[YY1], 0, 0), col = colZ1)
segments(x0 = 0, y0 = 0, x1 = max(x[xx]), y1 = 0, col = "white")
segments(x0 = 0, y0 = 0, x1 = 0, y1 = 1, col = "white")
seg.cat <- function(inter, cat, col.unf, col.exp, col.stab) {
if (inter[3] == 0) {
my.col = 0
}
if (inter[3] == 1) {
my.col = col.unf
}
if (inter[3] == 2) {
my.col = col.stab
}
if (inter[3] == -1) {
my.col = col.exp
}
segments(x0 = inter[1], y0 = -0.01, y1 = -0.01, x1 = inter[2],
col = my.col, lwd = 4, lty = 2)
}
cat <- ecospat.niche.dyn.index(z1, z2, intersection = quant)$dyn
inter <- cbind(z1$x[-length(z1$x)], z1$x[-1], cat[-1])
apply(inter, 1, seg.cat, col.unf = "#00FF0050",
col.exp = "#FF000050", col.stab = "#0000FF50")
}
if (!is.null(z1$y)) {
z <- t(as.matrix(z1$w + 2 * z2$w))[, nrow(as.matrix(z1$z.uncor)):1]
z1$Z <- t(as.matrix(z1$Z))[, nrow(as.matrix(z1$Z)):1]
z2$Z <- t(as.matrix(z2$Z))[, nrow(as.matrix(z2$Z)):1]
if (interest == 1) {
col1 = colorRampPalette(c("white", "red"),bias=1,alpha=0.8)(n=100)
col2 = colorRampPalette(c("white", "green"),bias=1)(n=100)
image(x = z1$x, y = z1$y, z = t(as.matrix(z1$z.uncor))[,
nrow(as.matrix(z1$z.uncor)):1], col = col1,
zlim = c(1e-05, cellStats(z1$z.uncor, "max")),
xlab = name.axis1, ylab = name.axis2)
image(x = z2$x, y = z2$y, z = t(as.matrix(z2$z.uncor))[,
nrow(as.matrix(z2$z.uncor)):1], col = col2,
zlim = c(1e-05, cellStats(z2$z.uncor, "max")), add = TRUE)
}
title(title)
contour(x = z1$x, y = z1$y, z1$Z, add = TRUE, levels = quantile(z1$Z[z1$Z >
0], c(0, quant)), drawlabels = FALSE, lty = c(1,
2), col = colZ1,lwd=3)
contour(x = z2$x, y = z2$y, z2$Z, add = TRUE, levels = quantile(z2$Z[z2$Z >
0], c(0, quant)), drawlabels = FALSE, lty = c(1,
2), col = colZ2,lwd=3)
contour(x = z2$x, y = z2$y, z, add = TRUE,col="gray65",lwd=1.5,drawlabels = FALSE)
}
}
arow <- function (sp1, sp2, clim1, clim2, col = "##FFFFFF00")
{
if (ncol(as.matrix(sp1)) == 2) {
arrows(median(sp1[, 1]), median(sp1[, 2]), median(sp2[,
1]), median(sp2[, 2]), col = "gray27", lwd = 4,
length = 0.1)
arrows(median(clim1[, 1]), median(clim1[, 2]), median(clim2[,
1]), median(clim2[, 2]), lty = "11", col = col,
lwd = 2, length = 0.1)
}
else {
arrows(median(sp1), 0.025, median(sp2), 0.025, col = "black",
lwd = 4, length = 0.1)
arrows(median(clim1), -0.025, median(clim2), -0.025,
lty = "11", col = col, lwd = 2, length = 0.1)
}
}